Preamble

This notebook deals with the comparison between inventories and fluxes of nitrate on one hand, and different aspects of primary production on the other, such as net primary production, new production, and export fluxes.

All figures exported from this notebook are prefixes by FIGURE_PP_.

In [1]:
%load_ext autoreload
%autoreload 2
%run imports.py

dims = dict(
    FN=hv.Dimension('FN', label='Nitrate flux', unit='mmol N m⁻² d⁻¹', range=(.005, 10)),
    pp_no3_eq=hv.Dimension('pp_no3_eq', label='NPP', unit='mmol N m⁻² d⁻¹'),
    flux=hv.Dimension('flux', label='N flux', unit='mmol N m⁻² d⁻¹', range=(.01, 10)),
)

Net primary production

First, we visualize annual net primary production inferred using ocean colour [@arrigo2015continued] as an average summer value during the last two decades. This map can then be compared with observed upward fluxes of nitrate during winter. It will not come as a surprise that NPP several times larger than the upward nitrate fluxes as it includes a considerable fraction regenerated production.

FIGURE_PP_map_NPP-FN

Annual average maps

Data were acquired as binary files through personal contact with G. van Dijken and K. Arrigo and converted to the netcdf format.

In [2]:
# agg = 'mean'
agg = 'summer_05-09_mean'
# agg = 'integral'
ds = xr.open_mfdataset(f'/Users/doppler/database/vandijken_arrigo_pp/annual/*{agg}.nc',
                      concat_dim='date').compute()
ds = ds.mean(dim='date')
df = ds.to_dataframe()
df['pp_no3_eq'] = df.pp / 12 * 16 / 106
# df.pp /= 1e3

/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/xarray/backends/api.py:783: FutureWarning: In xarray version 0.13 `auto_combine` will be deprecated.
  coords=coords)
/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/xarray/backends/api.py:783: FutureWarning: Also `open_mfdataset` will no longer accept a `concat_dim` argument.
To get equivalent behaviour from now on please use the new
`combine_nested` function instead (or the `combine='nested'` option to
`open_mfdataset`).The datasets supplied do not have global dimension coordinates. In
future, to continue concatenating without supplying dimension
coordinates, please use the new `combine_nested` function (or the
`combine='nested'` option to open_mfdataset.
  coords=coords)
/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/xarray/core/nanops.py:160: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)
In [3]:
def logmean(x):
    return np.nanmedian(np.log10(x))

options = (
    opts.HexTiles(aggregator=logmean, projection=ccrs.NorthPolarStereo(), gridsize=20,
                  tools=['hover'], colorbar_opts=dict(major_label_text_font_size='12pt'),
                  colorbar=True, width=570, height=510, cmap=default_cmap,
                  hooks=[logcolor_ticks([.4, 1., 3, 10, 30, 60])]),
    opts.Feature(scale='50m'),
)

pp = gv.HexTiles(df, vdims='pp_no3_eq')
l = pp * land
l.opts(*options)
hv.save(l.opts(toolbar=None, clone=True), #.redim.range(pp_no3_eq=(0.01, 2.5)),
              f'../nb_fig/FIGURE_PP_map_arrigo_vandijken_{agg}', fmt='png')
l

/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/numpy/lib/function_base.py:3405: RuntimeWarning: All-NaN slice encountered
  r = func(a, **kwargs)
/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/numpy/lib/function_base.py:3405: RuntimeWarning: All-NaN slice encountered
  r = func(a, **kwargs)
Out[3]:

Pan-Arctic patterns of annual net primary production show a very productive Atlantic sector, somewhat lower values in the Beaufort Sea, and a very oligotrophic central basin, similar to winter fluxes of nitrate. The main difference difference is that net primary production is at least an order of magnitude larger than any of the fluxes, so it seems hard to get more out of this than qualitative agreement at best.

The reason is of course that "net primary production" includes regenerated production, meaning the same nitrogen is just recycled over and over again, and can mostly do without input of new nitrogen.

In [4]:
display(HTML('../nb_fig/FIGURE_FN-COMP_map_winter.html'))

Pan-Arctic comparison of winter nitrate fluxes with annual new and export production

First, we need to compile pan-Arctic estimates of annual new

New production

Sakshaug's 2004 article in "The Organic Carbon Cycle in the Arctic Ocean" by Stein & MacDonald (eds.) has a good compilation. We use the following values:

In [5]:
gC_per_yr_to_mmolN_per_day = 1. /12 /106*16 * 1e3/365

data = {
    # Area: new production estimate in [gC / yr]
    'Barents Sea': [8., 100],
    'Chukchi Sea': [5, 160],
    'Beaufort Sea': [7, 17],
    'Baffin Bay': [25, 50],
    'Central basin': [.5]
}

df = pd.DataFrame(zip(*[(k, v) for k, vs in data.items() for v in vs]), index=['Area', 'newp']).transpose()
df.newp = gC_per_yr_to_mmolN_per_day * df.newp.astype(float)

df = df.groupby('Area').mean().reset_index()

df_newp = df.copy()

Export flux of particulate organic carbon

The basis of our brief discussion of POC export fluxes are the values compiled by I. Wiedmann's PhD thesis (2015, Fig. 2), supplemented with a number of additional references to increase data coverage in selected regions.

Ingrid Wiedmann, Fig. 2 PhD thesis

These values are all valid for a depth of 200 m.

Because the raw data to produce the graph was not available at the time of writing, we did some screen grabbing using a custom python script. This script is saved as grab_data.py, alongside the original and the reproduced figure.

In [6]:
# extracted raw coordinates relative to the png
xy = pd.read_csv('../data/wiedmann2015_export_fluxes/coords.csv', index_col=[0])[['x', 'y']]

# x0, ...: figure coordinates
# x0d, ...: data coordinates. y: logscale!
x0, x1, y0, y1 = xy.loc[0, 'x'], xy.loc[1, 'x'], xy.loc[0, 'y'], xy.loc[2, 'y']
x0d, x1d, y0d, y1d = 1., 12., 0., 3.

def convertx(c):
    return round((c-x0)/(x1-x0) * (x1d-x0d) + x0d)

def converty(c):
    return (c-y0)/(y1-y0) * (y1d-y0d) + y0d

df = dd.read_csv('../data/wiedmann2015_export_fluxes/wiedmann2015_fig2/*.csv').compute()

df['month'] = df.x.apply(convertx)
df['exp'] = 10**df.y.apply(converty)

df['exp_n_eq'] = df.exp /12 /106 *16
df = df.drop(columns=['x', 'y'])

# Sorting data

df['Season'] = 'Winter'
df.loc[df.month.between(4, 9), 'Season'] = 'Summer'

df['Area'] = ''

df.loc[df.ref.str.contains('barents'), 'Area'] = 'Barents Sea'
df.loc[df.ref.str.contains('nsvalbard'), 'Area'] = 'Barents Sea'
df.loc[df.ref.str.contains('beaufort-amundsen'), 'Area'] = 'Amundsen Gulf'
df.loc[df.ref.str.contains('baffin|nwp'), 'Area'] = 'Baffin Bay'
df.loc[df.ref.str.contains('kara'), 'Area'] = 'Kara'
df.loc[df.ref.str.contains('nwp'), 'Area'] = 'Northwater'
df.loc[df.ref.str.contains('greenland'), 'Area'] = 'Greenland Sea'
In [7]:
l = hv.Dataset(df, ['month', 'exp', 'ref']).to(hv.Points, groupby='ref').overlay().opts(
    opts.Points(legend_position='right', width=700, height=500, size=10, logy=True, padding=.1, jitter=.1)
)
hv.save(l.opts(toolbar=None), '../data/wiedmann2015_export_fluxes/Wiedmann-Figure2-recreated.png')

Some additional export estimates

These all pertain to the Canadian Basin.

In [8]:
#Honjo et al. 2010
flux = [0.05, 0.65] # gC m-2 yr-1

#Cai et al. 2010
flux.append(0.9) # gC m-2 yr-1

# convert to mmol N m-2 d-1
flux =  np.array(flux) / 12/106*16 *1e3/365

df = df.append(pd.DataFrame(dict(exp_n_eq=flux, Season='Summer', Area='Canadian Basin')), sort=False)

df_exp = df.copy()

Plot

In [9]:
options = opts.BoxWhisker(
    xrotation=45, xlabel='', logy=True, padding=.1, toolbar=None, yticks=[0.1, 1, 5, 10],
    height=500, width=500,
)

exp = hv.BoxWhisker(df_exp, ['Season', 'Area'], hv.Dimension('exp_n_eq', label='C export, Redfield N-equiv.', unit='mg N m-2 d-1'))
exp.opts(options)
Out[9]:

Nitrate fluxes

In [10]:
df_fn = pd.read_csv('../data/fn-compilation.csv')

fn = hv.BoxWhisker(df_fn, ['Season', 'Area'], hv.Dimension('FN', label='FN', unit='mg N m-2 d-1'))
fn.opts(options)

WARNING:param.BoxWhiskerPlot10326: Logarithmic axis range encountered value less than or equal to zero, please supply explicit lower-bound to override default of 0.055.
Out[10]:

Merge compilations of export production, new production, and nitrate fluxes

In [11]:
df_fn = df_fn.replace({'Perennial':'Winter', 'Barents Sea, AABC': 'Barents Sea',
                       'N Svalbard/Fram Strait': 'Barents Sea', 'Canada Basin': 'Central basin', 
                       'Makarov Basin': 'Central basin', 'Amundsen Gulf': 'Beaufort Sea'})

df_fn = df_fn.groupby(['Season', 'Area']).mean()['FN'].loc['Winter'].reset_index()

df_exp = df_exp.replace({'Amundsen Gulf': 'Beaufort Sea', 'Canadian Basin': 'Central basin'})

df_exp = df_exp.groupby(['Season', 'Area']).mean()['exp_n_eq'].loc['Summer'].reset_index()

df = df_fn.merge(df_exp, how='outer').merge(df_newp, how='outer').groupby('Area').mean()

df = df.stack().rename_axis(['Area', 'kind']).rename('flux').to_frame().reset_index()

df = df.replace({'exp_n_eq': 'Vertical export', 'FN': 'NO₃ flux', 'newp': 'New Production'})

df.head()
Out[11]:
Area kind flux
0 Vertical export 0.294907
1 Amundsen Basin NO₃ flux 0.050000
2 Baffin Bay NO₃ flux 1.700000
3 Baffin Bay Vertical export 0.253219
4 Baffin Bay New Production 1.292324
# Enter uncertainties manually from each paper # Bourgault et al, 2011 df_exp.loc[('Amundsen Gulf', 'NO3 flux'), 'flux_std'] = 0.04 # unpublished df_exp.loc[('Baffin Bay', 'NO3 flux'), 'flux_std'] = 0.5 # Randelhoff et al., 2015 df_exp.loc[('Barents Sea', 'NO3 flux'), 'flux_std'] = 0.5
In [12]:
hv.renderer('bokeh').theme = theme

opts_bk = [opts.Bars(
    color=hv.Cycle(['#AA5F77', '#2DA3A2', '#D6AE4A']),
    backend='bokeh',
    xrotation=45, xlabel='', logy=True, yticks=[.01, 0.1, 1, 10],
    width=800, height=400, 
    fontsize={v: 14 for v in ['ylabel', 'xlabel', 'yticks', 'xticks']},
    line_width=3, yformatter='%g', show_legend=True,
    tools=['hover'],
    toolbar=None,
)]
opts_mpl = translate_options(opts_bk, bokeh2mpl)

l = hv.Bars(df, ['Area', 'kind'], 'flux').opts(*opts_bk, *opts_mpl).redim(**dims)
l = l[['Baffin Bay', 'Barents Sea', 'Beaufort Sea', 'Central basin']]
fname = '../nb_fig/FIGURE_PP_FN_NEWP_EXP_regional'
hv.save(l, fname, fmt='png')
hv.save(l, fname, fmt='html')
save_bokeh_svg(l, fname+'.svg')
l
Out[12]:

Single case studies

First, manually enter the data from the studies.

In [34]:
## Randelhoff et al 2016 JGR

df = pd.DataFrame(dict(
    FN=[1.2, 0.6, 0.3, 1.1, 0.3, 0.1], 
    newP=[2.6, 3.1, 8.4, 0.015, 0.018, 0.048],
    #station=['P1', 'P3', 'P4', 'P5', 'P6', 'P7'], 
    season=['spring', 'spring', 'spring', 'summer', 'summer', 'summer']
))

df = df.set_index('season').stack().rename('flux')

df = df.rename_axis(['season', 'type']).reset_index()

df = df.replace(dict(spring='Spring', summer='Summer', newP='New prod.', FN='NO3 flux')).assign(ref='aRandelhoff et al.')

## Nishino et al 2018

df = df.append(
    pd.DataFrame({'New prod.': 0.24, 'NO3 flux': .19,}, index=[0]).stack()
    .rename('flux').rename_axis(('','type')).reset_index(1)
    .assign(season='Summer', ref='bNishino et al.'),
    sort=False, ignore_index=True
)
In [37]:
options = (
    opts.BoxWhisker(logy=True, xlabel='', xrotation=45, yformatter='%g',
                    width=600, height=330, box_fill_color='grey',
                    padding=.1, tools=['hover']),
)

l = hv.BoxWhisker(
    df, 
    [hv.Dimension('ref', value_format=lambda s: s[1:]), 'season', 'type'], 
    hv.Dimension('flux', label='N flux', range=(.01, 10), unit=u'mmol N m\u207B\u00B2 d\u207B\u00B9')
)
l = l.opts(*options)
fname = '../nb_fig/FIGURE_PP_incubations_vs_flux'
hv.save(l.opts(toolbar=None, clone=True), fname, fmt='png')
hv.save(l, fname, fmt='html')
l
Out[37]:

Paulsen et al., 2019

First, enter the values given in their Table 1.

In [32]:
df = pd.DataFrame(dict(
    Month=[1, 3, 5, 8, 11],
    DON=   [84., 51, 35, 109, 81],
    DONerr=[9, 32, 22, 36, 10],
    PON=   [4., 8, 46, 30, 3],
    PONerr=[1, 5, 41, 17, 1]
))

Then define an NdOverlay of Spread elements:

In [33]:
griddict = {pool: hv.Spread(df, 'Month', [pool, pool+'err'])
            for pool in ['PON', 'DON']}

options = [
    opts.Curve(padding=.1, line_width=3, width=500,
               xticks=[(1, 'Jan'), (3, 'Mar'), (5, 'May'), (8, 'Aug'), (11, 'Nov')],
               toolbar=None),
    opts.NdOverlay(legend_position='top_left'),
]

err = hv.HoloMap(griddict)
conc = err.map(hv.Curve, hv.Spread)

l = err*conc
l = l.overlay()
l.opts(*options)
for c, g in zip(['k', 'g'], ['PON', 'DON']):
    l[g].opts(opts.Curve(color=c), opts.Spread(color=c))

l = l.redim(DON=hv.Dimension('DON', label='DON, PON', unit='µM'))

fname = '../nb_fig/FIGURE_PP_chart_paulsen'
hv.save(l, fname, fmt='png')
hv.save(l, fname, fmt='html')
l
Out[33]: